Degradation statistics

baRulho:quantifying habitat-induced degradation of (animal) acoustic signals

Published

October 2, 2024

Source code, data and annotation protocol found at https://github.com/maRce10/barulho_paper

1 Purposes

Estimate degradation metrics on re-recorded signals from playback experiment at Bosque de Tlalpan, Mexico City, 2019
Run statistical analyses

Load packages

Code

# install/ load packages
sketchy::load_packages(packages = c("remotes", "kableExtra", "knitr",
    "formatR", "rprojroot", "maRce10/warbleR", "ggplot2", "tidyr",
    "viridis", "corrplot", "brms", "ggdist", "cowplot", "cmdstanr",
    "maRce10/brmsish"))

source("~/Dropbox/R_package_testing/brmsish/R/extended_summary.R")
source("~/Dropbox/R_package_testing/brmsish/R/helpers.R")

Code

degrad_df <- read.csv("./data/processed/tlalpan_degradation_metrics_v01.csv")

degrad_params <- c("blur.ratio", "spectrum.blur.ratio", "envelope.correlation",
    "excess.attenuation", "signal.to.noise.ratio", "cross.correlation",
    "tail.to.signal.ratio", "tail.to.noise.ratio", "spectrum.correlation")

comp.cases <- complete.cases(degrad_df[, names(degrad_df) %in% degrad_params])

pca <- prcomp(degrad_df[comp.cases, names(degrad_df) %in% degrad_params],
    scale. = TRUE)

# add to data
degrad_df$pc1 <- NA
degrad_df$pc1[comp.cases] <- pca$x[, 1]
degrad_df$pc1.1m.rate <- degrad_df$distance
# plot rotation values by PC
pca_rot <- as.data.frame(pca$rotation[, 1:4])

pca_rot_stck <- stack(pca_rot)

pca_rot_stck$variable <- rownames(pca_rot)
pca_rot_stck$Sign <- ifelse(pca_rot_stck$values > 0, "Positive", "Negative")
pca_rot_stck$rotation <- abs(pca_rot_stck$values)

ggplot(pca_rot_stck, aes(x = variable, y = rotation, fill = Sign)) +
    geom_col() + coord_flip() + scale_fill_viridis_d(alpha = 0.7,
    begin = 0.3, end = 0.8) + facet_wrap(~ind)

1.1 Correlation among metrics

Raw metrics:

Code

cormat <- cor(degrad_df[, degrad_params], use = "pairwise.complete.obs")

rownames(cormat) <- colnames(cormat) <- degrad_params

cols_corr <- colorRampPalette(c("white", "white", viridis(4, direction = -1)))(10)

cp <- corrplot.mixed(cormat, tl.cex = 0.7, upper.col = cols_corr,
    lower.col = cols_corr, order = "hclust", lower = "number", upper = "ellipse",
    tl.col = "black")

Code

# sort parameters as in clusters for cross correlation
degrad_params <- degrad_params[match(rownames(cp$corr), degrad_params)]

1.2 Data description

20 frequencies
3 locations
11520 test sounds
160 sound treatment combinations
Sample sizes per location, transect and signal type

Code

agg <- aggregate(cbind(sound.id, treatment.replicates) ~ location +
    habitat.structure + distance, degrad_df, function(x) length(unique(x)))

agg$replicates <- agg$sound.id/agg$treatment.replicates

agg

location	habitat.structure	distance	sound.id	treatment.replicates	replicates
1	closed	10	480	160	3
2	closed	10	480	160	3
3	closed	10	480	160	3
1	open	10	480	160	3
2	open	10	480	160	3
3	open	10	480	160	3
1	closed	30	480	160	3
2	closed	30	480	160	3
3	closed	30	480	160	3
1	open	30	480	160	3
2	open	30	480	160	3
3	open	30	480	160	3
1	closed	65	480	160	3
2	closed	65	480	160	3
3	closed	65	480	160	3
1	open	65	480	160	3
2	open	65	480	160	3
3	open	65	480	160	3
1	closed	100	480	160	3
2	closed	100	480	160	3
3	closed	100	480	160	3
1	open	100	480	160	3
2	open	100	480	160	3
3	open	100	480	160	3

1.3 Statistical analysis

Code

iter <- 10000
chains <- 4
priors <- c(prior(normal(0, 6), class = "b"), prior(cauchy(0, 4),
    class = "sd"))


# set frequency to mean-centered
degrad_df$frequency <- scale(degrad_df$frequency, scale = FALSE)

# set base level for factors
degrad_df$habitat.structure <- factor(degrad_df$habitat.structure,
    levels = c("open", "closed"))
degrad_df$frequency.modulation <- factor(degrad_df$frequency.modulation,
    levels = c("no_fm", "fm"))
degrad_df$amplitude.modulation <- factor(degrad_df$amplitude.modulation,
    levels = c("no_am", "am"))
degrad_df$duration <- factor(degrad_df$duration, levels = c("short",
    "long"))
degrad_df$location <- as.factor(degrad_df$location)
# degrad_df$ord.distance <- ordered(degrad_df$distance)

set.seed(123)

cmdstanr::set_cmdstan_path("~/Documentos/cmdstan/")

# to run within-chain parallelization
mod_pc1 <- brm(formula = pc1 ~ frequency * habitat.structure + frequency.modulation *
    habitat.structure + amplitude.modulation * habitat.structure +
    duration * habitat.structure + mo(distance) + (1 | location) +
    (1 | treatment.replicates), data = degrad_df, prior = priors,
    iter = iter, chains = chains, cores = chains, control = list(adapt_delta = 0.99,
        max_treedepth = 15), file = "./data/processed/regression_model_pc1.RDS",
    file_refit = "always", backend = "cmdstanr", threads = threading(10))

mod_blurratio <- brm(formula = blur.ratio ~ frequency * habitat.structure +
    frequency.modulation * habitat.structure + amplitude.modulation *
    habitat.structure + duration * habitat.structure + mo(distance) +
    (1 | location) + (1 | treatment.replicates), data = degrad_df,
    prior = priors, iter = iter, chains = chains, cores = chains,
    control = list(adapt_delta = 0.99, max_treedepth = 15), file = "./data/processed/regression_model_blur_ratio.RDS",
    file_refit = "always", backend = "cmdstanr", threads = threading(10))

Code

my.viridis <- function(...) viridis(alpha = 0.5, begin = 0.3, end = 0.7,
    ...)


effects <- c("habitat structure", "frequency modulation", "amplitude modulation",
    "frequency", "duration", "frequency:habitat structure", "habitat structure:duration",
    "habitat structure:amplitude modulation", "habitat structure:frequency modulation")
extended_summary(read.file = "./data/processed/regression_model_pc1.RDS",
    n.posterior = 1000, fill = "orange3", trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    )

1.4 regression_model_pc1

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 6) Intercept-student_t(3, -0.3, 2.5) sd-cauchy(0, 4) sigma-student_t(3, 0, 2.5) simo-dirichlet(1)	pc1 ~ frequency * habitat.structure + frequency.modulation * habitat.structure + amplitude.modulation * habitat.structure + duration * habitat.structure + mo(distance) + (1 \| location) + (1 \| treatment.replicates)	10000	4	1	5000	42 (0.0021%)	0	1795.606	3839.76	401962879

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
frequency	0.157	0.127	0.187	1.001	2607.552	4448.647
habitat structure	1.086	0.990	1.181	1	10835.391	12469.973
frequency modulation	0.050	-0.129	0.227	1.002	1795.606	4018.020
amplitude modulation	0.268	0.090	0.442	1.002	2014.067	3839.760
duration	-0.141	-0.320	0.033	1.002	2168.261	4800.424
frequency:habitat structure	0.047	0.031	0.063	1	19362.470	15343.730
habitat structure:frequency modulation	-0.087	-0.180	0.005	1	15136.276	15435.642
habitat structure:amplitude modulation	-0.007	-0.100	0.085	1	15877.418	14414.394
habitat structure:duration	-0.003	-0.097	0.092	1	15285.902	13945.763

Code

extended_summary(read.file = "./data/processed/regression_model_blur_ratio.RDS",
    n.posterior = 1000, fill = "orange3", trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    effects = effects)

1.5 regression_model_blur_ratio

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 6) Intercept-student_t(3, 0.1, 2.5) sd-cauchy(0, 4) sigma-student_t(3, 0, 2.5) simo-dirichlet(1)	blur.ratio ~ frequency * habitat.structure + frequency.modulation * habitat.structure + amplitude.modulation * habitat.structure + duration * habitat.structure + mo(distance) + (1 \| location) + (1 \| treatment.replicates)	10000	4	1	5000	100 (0.005%)	0	1905.433	4789.244	299027454

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
habitat structure	0.029	0.026	0.033	1	18137.687	14919.455
frequency modulation	0.065	0.058	0.072	1.002	2043.117	4835.342
amplitude modulation	0.023	0.016	0.030	1.002	1905.433	4789.244
frequency	0.001	0.000	0.002	1.001	2413.543	5026.225
duration	0.000	-0.006	0.008	1.003	2110.331	4830.679
frequency:habitat structure	0.003	0.002	0.003	1	20080.498	12256.359
habitat structure:duration	0.005	0.001	0.009	1	25094.631	15713.549
habitat structure:amplitude modulation	0.002	-0.001	0.006	1	25475.014	14808.890
habitat structure:frequency modulation	-0.020	-0.024	-0.017	1	27043.019	14561.904

Code

extended_summary(read.file = "./data/processed/regression_model_excess_attenuation.RDS",
    n.posterior = 1000, fill = "orange3", trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    effects = effects)

1.6 regression_model_excess_attenuation

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 6) Intercept-student_t(3, 81.7, 31.2) sd-cauchy(0, 4) sigma-student_t(3, 0, 31.2) simo-dirichlet(1)	excess.attenuation ~ frequency + frequency.modulation + amplitude.modulation + duration + habitat.structure + mo(distance_f) + (1 \| location) + (1 \| treatment.replicates)	20000	4	1	10000	5 (0.00012%)	0	3308.454	5800.994	1741877593

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
habitat structure	7.486	7.072	7.898	1	42076.716	29423.336
frequency modulation	-3.427	-4.985	-1.869	1.002	3308.454	7139.187
amplitude modulation	-3.295	-4.841	-1.686	1.001	3470.592	5800.994
frequency	13.703	12.116	15.299	1.001	3473.790	7645.888
duration	0.792	-0.824	2.390	1.001	3656.833	6971.910
NA	NA	NA	NA	NA	NA	NA
NA.1	NA	NA	NA	NA	NA	NA
NA.2	NA	NA	NA	NA	NA	NA
NA.3	NA	NA	NA	NA	NA	NA

Code

# extended_summary(read.file =
# './data/processed/regression_model_signal_to_noise_ratio.RDS',
# n.posterior = 1000, fill = 'orange3', trace.palette =
# my.viridis, remove.intercepts = TRUE, highlight = TRUE,
# gsub.pattern = c('b_', 'habitat.structureclosed',
# 'amplitude.modulationam', 'frequency.modulationfm',
# 'durationlong'), gsub.replacement = c('', 'habitat structure',
# 'amplitude modulation', 'frequency modulation', 'duration'),
# effects = effects)

extended_summary(read.file = "./data/processed/regression_model_tail_to_signal_ratio.RDS",
    n.posterior = 1000, fill = "orange3", trace.palette = my.viridis,
    remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_",
        "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm",
        "durationlong"), gsub.replacement = c("", "habitat structure",
        "amplitude modulation", "frequency modulation", "duration"),
    effects = effects)

1.7 regression_model_tail_to_signal_ratio

	priors	formula	iterations	chains	thinning	warmup	diverg_transitions	rhats > 1.05	min_bulk_ESS	min_tail_ESS	seed
1	b-normal(0, 6) Intercept-student_t(3, -6, 8.6) sd-cauchy(0, 4) sigma-student_t(3, 0, 8.6) simo-dirichlet(1)	tail.to.signal.ratio ~ frequency + frequency.modulation + amplitude.modulation + duration + habitat.structure + mo(distance_f) + (1 \| location) + (1 \| treatment.replicates)	20000	4	1	10000	26 (0.00065%)	0	3803.126	7149.158	1695545447

	Estimate	l-95% CI	u-95% CI	Rhat	Bulk_ESS	Tail_ESS
habitat structure	3.270	3.099	3.439	1	44650.292	29637.846
frequency modulation	-0.692	-1.297	-0.085	1.001	3901.160	7710.797
amplitude modulation	0.659	0.049	1.257	1	3803.126	8310.054
frequency	1.622	1.025	2.211	1	4186.279	8273.165
duration	-0.400	-1.006	0.202	1.002	3961.338	7149.158
NA	NA	NA	NA	NA	NA	NA
NA.1	NA	NA	NA	NA	NA	NA
NA.2	NA	NA	NA	NA	NA	NA
NA.3	NA	NA	NA	NA	NA	NA

2 Takeaways

Session information

R version 4.4.1 (2024-06-14)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.4 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=es_CR.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=es_CR.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=es_CR.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=es_CR.UTF-8 LC_IDENTIFICATION=C       

time zone: America/Costa_Rica
tzcode source: system (glibc)

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] brmsish_1.0.0       cmdstanr_0.8.1.9000 cowplot_1.1.3      
 [4] ggdist_3.3.2        brms_2.21.0         Rcpp_1.0.13        
 [7] corrplot_0.92       viridis_0.6.5       viridisLite_0.4.2  
[10] tidyr_1.3.1         ggplot2_3.5.1       warbleR_1.1.32     
[13] NatureSounds_1.0.4  seewave_2.2.3       tuneR_1.4.7        
[16] rprojroot_2.0.4     formatR_1.14        knitr_1.48         
[19] kableExtra_1.4.0    remotes_2.5.0      

loaded via a namespace (and not attached):
  [1] bitops_1.0-8         pbapply_1.7-2        gridExtra_2.3       
  [4] inline_0.3.19        testthat_3.2.1.1     sandwich_3.1-0      
  [7] rlang_1.1.4          magrittr_2.0.3       multcomp_1.4-25     
 [10] matrixStats_1.3.0    compiler_4.4.1       loo_2.8.0           
 [13] reshape2_1.4.4       systemfonts_1.1.0    vctrs_0.6.5         
 [16] stringr_1.5.1        arrayhelpers_1.1-0   pkgconfig_2.0.3     
 [19] crayon_1.5.3         fastmap_1.2.0        backports_1.5.0     
 [22] labeling_0.4.3       fontawesome_0.5.2    utf8_1.2.4          
 [25] rmarkdown_2.28       ps_1.7.7             purrr_1.0.2         
 [28] xfun_0.47            jsonlite_1.8.8       highr_0.11          
 [31] tidybayes_3.0.6      uuid_1.2-1           parallel_4.4.1      
 [34] R6_2.5.1             stringi_1.8.4        StanHeaders_2.32.10 
 [37] brio_1.1.5           estimability_1.5.1   rstan_2.32.6        
 [40] zoo_1.8-12           bayesplot_1.11.1     Matrix_1.7-0        
 [43] splines_4.4.1        tidyselect_1.2.1     rstudioapi_0.16.0   
 [46] abind_1.4-5          yaml_2.3.10          dtw_1.23-1          
 [49] codetools_0.2-20     processx_3.8.4       curl_5.2.2          
 [52] pkgbuild_1.4.4       plyr_1.8.9           lattice_0.22-6      
 [55] tibble_3.2.1         withr_3.0.1          bridgesampling_1.1-2
 [58] posterior_1.6.0      coda_0.19-4.1        evaluate_0.24.0     
 [61] signal_1.8-1         survival_3.7-0       sketchy_1.0.3       
 [64] proxy_0.4-27         RcppParallel_5.1.8   xml2_1.3.6          
 [67] pillar_1.9.0         tensorA_0.36.2.1     packrat_0.9.2       
 [70] checkmate_2.3.2      stats4_4.4.1         distributional_0.4.0
 [73] generics_0.1.3       RCurl_1.98-1.16      rstantools_2.4.0    
 [76] munsell_0.5.1        scales_1.3.0         xtable_1.8-4        
 [79] glue_1.7.0           emmeans_1.10.3       tools_4.4.1         
 [82] xaringanExtra_0.8.0  mvtnorm_1.2-5        grid_4.4.1          
 [85] ape_5.8              QuickJSR_1.3.1       colorspace_2.1-1    
 [88] nlme_3.1-165         cli_3.6.3            svUnit_1.0.6        
 [91] fansi_1.0.6          svglite_2.1.3        Brobdingnag_1.2-9   
 [94] dplyr_1.1.4          V8_4.4.2             gtable_0.3.5        
 [97] fftw_1.0-8           digest_0.6.37        TH.data_1.1-2       
[100] farver_2.1.2         rjson_0.2.22         htmlwidgets_1.6.4   
[103] htmltools_0.5.8.1    lifecycle_1.0.4      MASS_7.3-61

--- title: Degradation statistics subtitle: baRulho:quantifying habitat-induced degradation of (animal) acoustic signals # author: Marcelo Araya-Salas date: "`r Sys.Date()`" toc: true toc-depth: 2 toc-location: left number-sections: true highlight-style: pygments format: html: df-print: kable code-fold: true code-tools: true code-copy: true embed-resources: true editor_options: chunk_output_type: console ---   ```{r set root directory, echo = FALSE} knitr::opts_knit$set(root.dir = "..") ``` ```{r add link to github repo, echo = FALSE, results='asis'} # print link to github repo if any if (file.exists("./.git/config")){ config <- readLines("./.git/config") url <- grep("url", config, value = TRUE) url <- gsub("\\turl = |.git$", "", url) cat("\nSource code, data and annotation protocol found at [", url, "](", url, ")", sep = "") } ``` ```{r setup style, echo = FALSE, message = FALSE, warning=FALSE} # options to customize chunk outputs knitr::opts_chunk$set( tidy.opts = list(width.cutoff = 65), tidy = TRUE, message = FALSE, warning = FALSE ) ``` <div class="alert alert-info"> # Purposes - Estimate degradation metrics on re-recorded signals from playback experiment at Bosque de Tlalpan, Mexico City, 2019 - Run statistical analyses </div>   # Load packages {.unnumbered .unlisted} ```{r load packages} # install/ load packages sketchy::load_packages(packages = c("remotes", "kableExtra", "knitr", "formatR", "rprojroot", "maRce10/warbleR", "ggplot2", "tidyr", "viridis", "corrplot", "brms", "ggdist", "cowplot", "cmdstanr", "maRce10/brmsish")) source("~/Dropbox/R_package_testing/brmsish/R/extended_summary.R") source("~/Dropbox/R_package_testing/brmsish/R/helpers.R") ``` ```{r load packages and setup style, echo = FALSE, message = FALSE, warning=FALSE} # github packages must include user name ("user/package") # knitr is require for creating html/pdf/word reports # kableExtra is used to print pretty formatted tables # formatR is used for soft-wrapping code # klippy is used for adding a copy button to each code block pkgs <- c("remotes", "kableExtra", "knitr", "formatR", "rprojroot", "maRce10/warbleR", "ggplot2", "tidyr", "viridis", "corrplot", "brms", "ggdist", "cowplot", "cmdstanr", "maRce10/brmsish") # install/ load packages out <- lapply(pkgs, function(y) { # get pakage name pkg <- strsplit(y, "/")[[1]] pkg <- pkg[length(pkg)] # check if installed, if not then install if (!pkg %in% installed.packages()[,"Package"]) { if (grepl("/", y)) remotes::install_github(y, force = TRUE) else install.packages(y) } # load package a <- try(require(pkg, character.only = T), silent = T) if (!a) remove.packages(pkg) }) # set working directory as project directory or one directory above, rootdir <- try(rprojroot::find_rstudio_root_file(), silent = TRUE) if (is(rootdir, "try-error")) rootdir <- ".." opts_knit$set(root.dir = rootdir) # options to customize chunk outputs knitr::opts_chunk$set( class.source = "numberLines lineAnchors", # for code line numbers tidy.opts = list(width.cutoff = 65), tidy = TRUE, message = FALSE ) # this is a customized printing style data frames # screws up tibble function tibble <- function(x, ...) { x <- kbl(x, digits=4, align= 'c', row.names = FALSE) x <- kable_styling(x, position ="center", full_width = FALSE, bootstrap_options = c("striped", "hover", "condensed", "responsive")) asis_output(x) } registerS3method("knit_print", "data.frame", tibble) # to add copy button to code blocks htmltools::tagList( xaringanExtra::use_clipboard( button_text = "<i class=\"fa fa-clipboard\"></i>", success_text = "<i class=\"fa fa-check\" style=\"color: #90BE6D\"></i>", error_text = "<i class=\"fa fa-times-circle\" style=\"color: #F94144\"></i>" ), rmarkdown::html_dependency_font_awesome() ) theme_set(theme_classic(base_size = 17)) ``` ```{r read data degradation, eval = TRUE, echo = FALSE, warning=FALSE, message=FALSE, include=FALSE} degrad_df <- read.csv("./data/processed/barulho_degradation_metrics.csv", stringsAsFactors = FALSE) degrad_params <- c(bl.rt = "blur.ratio", sp.bl.rt = "spetrum.blur.ratio", env.cr = "envelope.correlation", EA ="excess.attenuation", SPL = "SPL", SNR = "signal.to.noise.ratio", SPCC = "cross.correlation", TSR = "tail.to.signal.ratio", TNR = "tail.to.noise.ratio", sp.cr = "spectrum.correlation", pc1 = "pc1") master_metadata <- read.csv("./data/raw/consolidated.master.sf_selection_table.csv") # keep only simulated data (exluding data from a student) master_metadata <- master_metadata[2:961, ] # infer original frequency possb_freqs <- seq(0.5, 10, length.out = 20) master_metadata$freq <- sapply(1:nrow(master_metadata), function(x) possb_freqs[which.min(abs((master_metadata$top.freq[x] + master_metadata$bottom.freq[x]) / 2 - possb_freqs))] ) freq_contour <- freq_ts(master_metadata, length.out = 20, parallel = 22, img = FALSE, path = "./data/raw/master_sound_file", pb = FALSE) freq_contour$mean.freq <- rowMeans(freq_contour[, 3:ncol(freq_contour)]) freq_contour$orig_sound_file <- master_metadata$orig.sound.file degrad_df$frequency <- sapply(1:nrow(degrad_df), function(x) freq_contour$mean.freq[freq_contour$orig_sound_file == degrad_df$sound.id[x]] ) degrad_df$sim.frequency <- sapply(1:nrow(degrad_df), function(x) master_metadata$freq[master_metadata$orig.sound.file == degrad_df$sound.id[x]] ) degrad_df <- separate(degrad_df, col = sound.id, into = c("Duration", "Freq_modulation", "Amp_modulation", "Harmonicity", "delete"), remove = FALSE, sep = "_") degrad_df$Duration <- ifelse(degrad_df$Duration == 0.1, "short", "long") degrad_df$`Freq_modulation` <- ifelse(degrad_df$`Freq_modulation` == "BB", "fm", "no_fm") degrad_df$Harmonicity <- gsub(".wav", "", degrad_df$Harmonicity) degrad_df$delete <- NULL degrad_df$treatment.replicates <- paste0("frq:", "_", degrad_df$sim.frequency, "dur:", degrad_df$Duration, "_", degrad_df$Amp_modulation, "_", degrad_df$Freq_modulation) degrad_df$Freq_modulation <- factor(degrad_df$Freq_modulation, levels = c("no_fm", "fm")) degrad_df$Amp_modulation[degrad_df$Amp_modulation == "no.am"] <- "no_am" degrad_df$Amp_modulation <- factor(degrad_df$Amp_modulation, levels = c("no_am", "am")) degrad_df$duration <- factor(degrad_df$Duration, levels = c("short", "long")) degrad_df <- degrad_df[degrad_df$Harmonicity != "harm", ] degrad_df <- degrad_df[degrad_df$distance > 1, ] degrad_df$frequency <- degrad_df$sim.frequency degrad_df$duration <- degrad_df$Duration degrad_df$amplitude.modulation <- degrad_df$Amp_modulation degrad_df$frequency.modulation <- degrad_df$Freq_modulation degrad_df$location <- degrad_df$Transecto # degrad_df$spectrum.blur.ratio <- degrad_df$spetral.blur.ratio degrad_df$habitat.structure <- ifelse(degrad_df$Vegetacion == "abierta", "open", "closed") degrad_df$sim.frequency <- degrad_df$signal.to.noise.ratio.t2 <- degrad_df$SPL <- degrad_df$template <- degrad_df$data.set <- degrad_df$Temperatura <- degrad_df$Harmonicity <- degrad_df$freq <- degrad_df$Humedad <- degrad_df$Vegetacion <- degrad_df$Duration <- degrad_df$orig.sound.file <- degrad_df$Freq_modulation <- degrad_df$Amp_modulation <- degrad_df$spetral.blur.ratio <- degrad_df$Transecto <- NULL names(degrad_df) degrad_df <- degrad_df[, c("sound.files", "selec", "start", "end", "bottom.freq", "top.freq", "sound.id", "habitat.structure", "distance", "reference", "frequency", "duration", "treatment.replicates", "location", "frequency.modulation", "amplitude.modulation","blur.ratio", "envelope.correlation", "excess.attenuation", "signal.to.noise.ratio", "cross.correlation", "spectrum.blur.ratio", "spectrum.correlation", "tail.to.noise.ratio", "tail.to.signal.ratio")] write.csv(degrad_df, "./data/processed/tlalpan_degradation_metrics_v01.csv", row.names = FALSE) ``` ```{r add PCA, out.width = "100%"} degrad_df <- read.csv("./data/processed/tlalpan_degradation_metrics_v01.csv") degrad_params <- c("blur.ratio", "spectrum.blur.ratio", "envelope.correlation", "excess.attenuation", "signal.to.noise.ratio", "cross.correlation", "tail.to.signal.ratio", "tail.to.noise.ratio", "spectrum.correlation") comp.cases <- complete.cases(degrad_df[,names(degrad_df) %in% degrad_params]) pca <- prcomp(degrad_df[comp.cases, names(degrad_df) %in% degrad_params], scale. = TRUE) # add to data degrad_df$pc1 <- NA degrad_df$pc1[comp.cases] <- pca$x[, 1] degrad_df$pc1.1m.rate <- degrad_df$distance # plot rotation values by PC pca_rot <- as.data.frame(pca$rotation[, 1:4]) pca_rot_stck <- stack(pca_rot) pca_rot_stck$variable <- rownames(pca_rot) pca_rot_stck$Sign <- ifelse(pca_rot_stck$values > 0, "Positive", "Negative") pca_rot_stck$rotation <- abs(pca_rot_stck$values) ggplot(pca_rot_stck, aes(x = variable, y = rotation, fill = Sign)) + geom_col() + coord_flip() + scale_fill_viridis_d(alpha = 0.7, begin = 0.3, end = 0.8) + facet_wrap(~ ind) ``` ## Correlation among metrics Raw metrics: ```{r correlation among raw metrics, out.width = "120%"} cormat <- cor(degrad_df[, degrad_params], use = "pairwise.complete.obs") rownames(cormat) <- colnames(cormat) <- degrad_params cols_corr <- colorRampPalette(c("white", "white", viridis(4, direction = -1)))(10) cp <- corrplot.mixed(cormat, tl.cex = 0.7, upper.col = cols_corr, lower.col = cols_corr, order = "hclust", lower = "number", upper = "ellipse", tl.col = "black") # sort parameters as in clusters for cross correlation degrad_params <- degrad_params[match(rownames(cp$corr), degrad_params)] ``` ## Data description - `r length(unique(degrad_df$frequency))` frequencies - `r length(unique(degrad_df$location))` locations - `r nrow(degrad_df)` test sounds - `r length(unique(degrad_df$treatment.replicates))` sound treatment combinations - Sample sizes per location, transect and signal type ```{r} agg <- aggregate(cbind(sound.id, treatment.replicates) ~ location + habitat.structure + distance, degrad_df, function(x) length(unique(x))) agg$replicates <- agg$sound.id / agg$treatment.replicates agg ``` ## Statistical analysis ```{r, eval =FALSE} iter <- 10000 chains <- 4 priors <- c(prior(normal(0, 6), class = "b"), prior(cauchy(0, 4), class = "sd")) # set frequency to mean-centered degrad_df$frequency <- scale(degrad_df$frequency, scale = FALSE) # set base level for factors degrad_df$habitat.structure <- factor(degrad_df$habitat.structure, levels = c("open", "closed")) degrad_df$frequency.modulation <- factor(degrad_df$frequency.modulation, levels = c("no_fm", "fm")) degrad_df$amplitude.modulation <- factor(degrad_df$amplitude.modulation, levels = c("no_am", "am")) degrad_df$duration <- factor(degrad_df$duration, levels = c("short", "long")) degrad_df$location <- as.factor(degrad_df$location) # degrad_df$ord.distance <- ordered(degrad_df$distance) set.seed(123) cmdstanr::set_cmdstan_path("~/Documentos/cmdstan/") # to run within-chain parallelization mod_pc1 <- brm( formula = pc1 ~ frequency * habitat.structure + frequency.modulation * habitat.structure + amplitude.modulation * habitat.structure + duration * habitat.structure + mo(distance) + (1 | location) + (1 | treatment.replicates), data = degrad_df, prior = priors, iter = iter, chains = chains, cores = chains, control = list(adapt_delta = 0.99, max_treedepth = 15), file = "./data/processed/regression_model_pc1.RDS", file_refit = "always", backend = "cmdstanr", threads = threading(10) ) mod_blurratio <- brm( formula = blur.ratio ~ frequency * habitat.structure + frequency.modulation * habitat.structure + amplitude.modulation * habitat.structure + duration * habitat.structure + mo(distance) + (1 | location) + (1 | treatment.replicates), data = degrad_df, prior = priors, iter = iter, chains = chains, cores = chains, control = list(adapt_delta = 0.99, max_treedepth = 15), file = "./data/processed/regression_model_blur_ratio.RDS", file_refit = "always", backend = "cmdstanr", threads = threading(10) ) ``` ```{r, eval = TRUE, results='asis', warning = FALSE, out.width="100%"} my.viridis <- function(...) viridis(alpha = 0.5, begin = 0.3, end = 0.7, ...) effects <- c("habitat structure", "frequency modulation", "amplitude modulation", "frequency", "duration", "frequency:habitat structure", "habitat structure:duration", "habitat structure:amplitude modulation", "habitat structure:frequency modulation") extended_summary(read.file = "./data/processed/regression_model_pc1.RDS", n.posterior = 1000, fill = "orange3", trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm", "durationlong"), gsub.replacement = c("", "habitat structure", "amplitude modulation", "frequency modulation", "duration"), ) extended_summary(read.file = "./data/processed/regression_model_blur_ratio.RDS", n.posterior = 1000, fill = "orange3", trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm", "durationlong"), gsub.replacement = c("", "habitat structure", "amplitude modulation", "frequency modulation", "duration"), effects = effects) extended_summary(read.file = "./data/processed/regression_model_excess_attenuation.RDS", n.posterior = 1000, fill = "orange3", trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm", "durationlong"), gsub.replacement = c("", "habitat structure", "amplitude modulation", "frequency modulation", "duration"), effects = effects) # extended_summary(read.file = "./data/processed/regression_model_signal_to_noise_ratio.RDS", # n.posterior = 1000, fill = "orange3", trace.palette = my.viridis, # remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_", # "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm", # "durationlong"), gsub.replacement = c("", "habitat structure", # "amplitude modulation", "frequency modulation", "duration"), effects = effects) extended_summary(read.file = "./data/processed/regression_model_tail_to_signal_ratio.RDS", n.posterior = 1000, fill = "orange3", trace.palette = my.viridis, remove.intercepts = TRUE, highlight = TRUE, gsub.pattern = c("b_", "habitat.structureclosed", "amplitude.modulationam", "frequency.modulationfm", "durationlong"), gsub.replacement = c("", "habitat structure", "amplitude modulation", "frequency modulation", "duration"), effects = effects) ```  <div class="alert alert-success">   # Takeaways </div>    ---    <font size="4">Session information</font> ```{r session info, echo=F} sessionInfo() ```